﻿using System;
using System.Collections.Generic;
using LatoolNet;

namespace NumericalSolution {
  class Diffusion1DFEM {
    private int m_numberOfElements;
    private int m_numberOfNodes;
    private int m_numberOfNodesInAnElement;
    private NumericalIntegration m_integration;
    private Element[] m_elements;
    private Matrix m_A;
    private Matrix m_B;
    private double[] m_initial_u;
    private double m_deltaTime;
    private double m_currentTime;
    private double m_endTime;
    private double m_deltaX;
    private double m_lambda;
    private double[] m_vec;

    public int NumberOfElements {
      get { return m_numberOfElements; }
      set {
        m_numberOfElements = value;
        m_numberOfNodes = value + 1;
      }
    }

    public int NumberOfNodes {
      get { return m_numberOfNodes; }
    }

    public int NumberOfNodesInAnElement {
      get { return m_numberOfNodesInAnElement; }
      set { m_numberOfNodesInAnElement = value; }
    }

    public double DeltaX {
      get { return m_deltaX; }
      set { m_deltaX = value; }
    }

    public double DeltaTime {
      get { return m_deltaTime; }
      set { m_deltaTime = value; }
    }

    public double CurrentTime {
      get { return m_currentTime; }
      set { m_currentTime = value; }
    }

    public double EndTime {
      get { return m_endTime; }
      set { m_endTime = value; }
    }

    public NumericalIntegration NumericalIntegration {
      get { return m_integration; }
      set { m_integration = value; }

    }
    public void CreateMesh() {

      m_elements = new Element[m_numberOfElements];
      int number = 0;
      for (int i = 0; i < m_numberOfElements; i++) {
        Element e = new Element(m_numberOfNodesInAnElement);
        m_elements[i] = e;
        e.Number = number++;
        e.Node[0].X = i * m_deltaX;
        e.Node[1].X = (i + 1) * m_deltaX;
      }
    }

    public void Solve() {

      m_A = new Matrix(m_numberOfNodes, m_numberOfNodes, 
                       MatrixType.DoubleSymmetric);
      m_B = new Matrix(m_numberOfNodes, m_numberOfNodes, 
                       MatrixType.DoubleSymmetric);

      for (int i = 0; i < m_numberOfElements; i++) {
        Matrix A_e = new Matrix(m_numberOfNodesInAnElement,
                                m_numberOfNodesInAnElement);
        Matrix B_e = new Matrix(m_numberOfNodesInAnElement,
                                m_numberOfNodesInAnElement);

        Element element = m_elements[i];
        AssembleLocalMatrix(A_e, B_e, element);
        AssembleGlobalMatrix(A_e, B_e, element);
      }

      //EigenvalueProblem p = new EigenvalueProblem();
      //p.GeneralizedSolve(m_A, m_B);

      //m_lambda = p.EigenValue[1, 0];

      //m_vec = new double[p.EigenVector.ColNum];
      //for (int i = 0; i < p.EigenVector.ColNum - 1; i++) {
      //  m_vec[i] = Math.Abs(p.EigenVector[i, 0]);
      //}

      double[] resultValues;
      Matrix resultVectors;
      Eigenproblem.Solve(m_A, m_B, out resultValues, out resultVectors);

      m_lambda = resultValues[1];

      m_vec = new double[resultVectors.RowNum];
      for (int i = 0; i < resultVectors.RowNum; i++) {
        m_vec[i] = Math.Abs(resultVectors[i, 0]);
      }

    }

    private void AssembleLocalMatrix(Matrix A_e, 
                                     Matrix B_e, 
                                     Element element) {

      for (int integrationIndex = 0;
           integrationIndex < m_integration.NumberOfIntegrationPoints;
           integrationIndex++) {

        Integrands(A_e, B_e, element,
                   m_integration.IntegrationPoint(integrationIndex),
                   m_integration.Weight(integrationIndex));
      }
    }

    private void Integrands(Matrix A_e, Matrix B_e, Element element,
                            double xi, double weight) {
      double h = element.Node[1].X - element.Node[0].X;
      double J = h / 2;

      for (int r = 0; r < m_numberOfNodesInAnElement; r++) {
        for (int s = 0; s < m_numberOfNodesInAnElement; s++) {
          A_e[r, s] += dPhi(r, xi) * (2 / h) *
                       dPhi(s, xi) * (2 / h) *
                       J * weight;
          B_e[r, s] += Phi(r, xi) * Phi(s, xi) * J * weight;
        }
      }
    }

    private void AssembleGlobalMatrix(Matrix A_e, Matrix B_e, Element e) {

      for (int r = 0; r < m_numberOfNodesInAnElement; r++) {
        for (int s = 0; s < m_numberOfNodesInAnElement; s++) {
          if (r > s) {
            //this condition is needed because symmetric elements are updated
            //at the same time. For example, if a[i,j] is updated then a[j,i]
            //is updated automatically by the LatoolNet library.
          } else {
            m_A[q(e.Number, r), q(e.Number, s)] += A_e[r, s];
            m_B[q(e.Number, r), q(e.Number, s)] += B_e[r, s];
          }
        }
      }
    }

    public double[] Initialize() {
      m_currentTime = 0.0;

      m_initial_u = new double[m_numberOfNodes];
      m_initial_u[0] = 0.0;
      for (int i = 1; i < m_numberOfNodes - 1; i++) {
        m_initial_u[i] = Math.Sin(Math.PI * i * m_deltaX);
      }
      m_initial_u[m_numberOfNodes - 1] = 0.0;

      return m_initial_u;
    }

    public double[] Next() {

      double[] current_u = new double[m_numberOfNodes];
      current_u[0] = 0.0;
      for (int i = 0; i < m_numberOfNodes; i++) {
        current_u[i] = m_initial_u[i] 
          * m_vec[i] * Math.Exp(- m_lambda * m_currentTime);
      }
      current_u[m_numberOfNodes - 1] = 0.0;

      return current_u;
     
    }

    private double Phi(int type, double xi) {
      if (type == 0) {
        return 0.5 * (1 - xi);
      } else if (type == 1) {
        return 0.5 * (1 + xi);
      }
      throw new ArgumentOutOfRangeException();
    }

    private double dPhi(int type, double xi) {
      if (type == 0) {
        return -0.5;
      } else if (type == 1) {
        return 0.5;
      }
      throw new ArgumentOutOfRangeException();
    }

    private int q(int elementNumber, int localNumber) {
      return elementNumber + localNumber;
    }
  }
}
